Velocity and energy distributions in microcanonical ensembles of hard spheres 
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In a microcanonical ensemble (constant NVE, hard reflecting walls) and in a molecular dynamics 
ensemble (constant NVEPG, periodic boundary conditions) with a number N of hard spheres in 
a d-dimensional volume V having a total energy E, a total momentum P, and an overall center of 
mass position G, the individual velocity components, velocity moduli, and energies follow trans- 
formed beta distributions with different arguments and shape parameters depending on d, TV, E, 
the boundary conditions, and possible symmetries in the initial conditions. This can be shown 
marginalizing the joint distribution of individual energies, which is a symmetric Dirichlet distri- 
bution. In the thermodynamic limit the beta distributions converge to gamma distributions with 
different arguments and shape or scale parameters, corresponding respectively to the Gaussian, i.e., 
Maxwell-Boltzmann, Maxwell, and Boltzmann or Boltzmann-Gibbs distribution. These analyti- 
cal results are in agreement with molecular dynamics and Monte Carlo simulations with different 
numbers of hard disks or spheres and hard reflecting walls or periodic boundary conditions. 

PACS numbers: 02.50.-r 02.50.Ng, 02.70. Uu, 05.10.-a, 05.10.Ln, 07.05.Tp 



I. INTRODUCTION 

The problem of the velocity distribution in a gas of 
hard spheres was discussed in a paper published by 
Maxwell in 1860 [![. Maxwell obtained the velocity dis- 
tribution by assuming independence of the three compo- 
nents of velocity and rotational invariance of the joint 
distribution. The only distribution satisfying the func- 
tional equation 



f v (xi,x 2 ,x 3 ) 



has factors of the form 



c 2 + x a) 

fv 1 ( x l)fv 2 (x2)fv 3 (x 3 ) 



f Va (x) = Aexp(-Bx 2 ), 



(1) 



(2) 



a = 1,2,3. This simple heuristic derivation can still 
be found in modern textbooks in statistical physics or 
physical chemistry [2j, but generalizations of Maxwell's 
method appeared earlier in the physical literature |3J . 

In 1867 Maxwell [|| became aware that Eq. should 
appear as a stationary solution for the dynamics of 
the gas and introduced a concept that later was called 
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Stofizahlansatz by Boltzmann. It led to a more detailed 
study of molecular collisions and to kinetic equations 
whose stationary solutions coincide with Maxwell's orig- 
inal distribution (see Refs. [5|-|9| for a modern mathe- 
matical approach to kinetic equations). This route was 
followed by Boltzmann, who obtained the velocity dis- 
tribution in a more general way in a series of papers 
written between 1868 and 1871 [I3-[ll. Based on the 
Stofizahlansatz, Boltzmann could prove that Maxwell's 
distribution is stationary. These results are summarized 
in Tolman's book [l~3| and in the first chapter of ter 
Haar's book [3~3 | . 

The two physicists were not working in isolation and 
were aware of their respective works. In his 1872 paper, 
Boltzmann often quotes Maxwell [ill . In 1873, Maxwell 
wrote to his correspondent Tait [16| : 

By the study of Boltzmann I have been un- 
able to understand him. He could not under- 
stand me on account of my shortness, and his 
length was and is an equal stumbling block to 
me. 

More details on the relationship between Maxwell and 
Boltzmann and on the influence of Maxwell on Boltz- 
mann's thought have been collected by Uffink [l7j . 

Tolman's analysis of classical binary collisions for hard 
spheres led to rate equations which can be interpreted as 
transition probabilities for a Markov chain after proper 
normalization. The interested reader can consult chap- 
ter V of Tolman's classic book [13J, in particular the 
discussion around Eq. (45.3) on page 129. The connec- 
tion with Markov chains was made explicit by Costan- 
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tini and Garibaldi [l8|, Ojj], who used a model due to 
Brillouin [20|. Before Costantini and Garibaldi, Penrose 
suggested that a Markovian hypothesis could justify the 
use of standard statistical mechanical tools [21j • Accord- 
ing to our interpretation of Penrose, due to the limits 
in human knowledge naturally leading to coarse grain- 
ing, systems of many interacting particles effectively be- 
have as Markov chains. Moreover, the possible number 
of states of such a chain is finite even if very large, there- 
fore only the theory of finite Markov chains is useful. 
Statistical equilibrium is reached when the system states 
obey the equilibrium distribution of the finite Markov 
chain; this equilibrium distribution exists, is unique and 
coincides with the stationary distribution if the chain is 
irreducible or ergodic. This point of view is also known as 
Markovianism. Indeed, in a recent paper on the Ehren- 
fest urn, we showed that, after appropriate coarse grain- 
ing, a Markov chain well approximates the behaviour of 
a realistic model for a fluid |22j |. 

Here we study the velocity distribution in a system of 
N smooth elastic hard spheres in d dimensions. Even if 
the evolution of the system is deterministic, we can con- 
sider the velocity components of each particle as random 
variables. We do not consider a finitary [23| version of 
the model by discretizing velocities, but keep them as real 
variables. Then a heuristic justification of Eq. © can be 
based on the central limit theorem (CLT). Here is the ar- 
gument. Following Maxwell's idea, one can consider the 
velocity components of each particle independent from 
each other. Further assuming that velocity jumps after 
collisions are independent and identically distributed ran- 
dom variables, one obtains for the velocity component a 
of a particle i at time t 

n(t) 

Via (t) = Vj a (0) + Av laJ , (3) 

i=i 

where n(t) is the number of collisions for that particle up 
to time t and Avi a j is the change in velocity at collision 
j. If the hypotheses stated above are valid, Eq. © de- 
fines a continuous-time random walk and the distribution 
fv ia (x, t) approaches a normal distribution for large t as 
a consequence of the CLT. Unfortunately this argument 
is only approximately true in the case of large systems 
and false for smaller systems. 

In Sec. |H] we obtain the theoretical probability den- 
sity functions of the individual energies, velocity moduli 
and velocity components, starting from the fundamental 
uniform distribution law in phase space. In Sec. Mil we 
present the molecular dynamics method used to simulate 
hard spheres. Interestingly, the same distributions can be 
reproduced by a simple Monte Carlo stochastic model in- 
troduced in Sec. IIVI The numerical results are presented 
in Sec. [V] together with some statistical goodness-of-fit 
tests. Indeed, it turns out that an equilibrium distri- 
bution of the velocity components seems to be reached 
already for N — 2 particles and without using any coarse 
graining. When N grows the equilibrium distribution ap- 



proaches the normal distribution, Eq. ([2]). A discussion 
and a summary follow in Sec. IVII 

II. THEORY 

We consider a fluid of N hard spheres in d dimensions 
with the same diameter a and mass m in a cuboidal 
box with sides L a , a — l,...,d. The positions r,, 
i = 1 , . . . , N are confined to a <i-dimensional box with 
volume V = ria=i -^ Q ' ^ e '' eacn position component ri a 
can vary in the interval [—L a /2, L a /2]. Elastic collisions 
transfer kinetic energy between the particles, while the 
total energy of the system, 

1 N 1 
E = - ^2 ■ Vj = -m v • v, (4) 

i=i 

does not change in time, i.e., it is a constant of the mo- 
tion. Therefore, the velocities v, are confined to the sur- 
face of a hypersphere given by the constraint that the 
total energy is E, i.e., each velocity component vi a can 
vary in the interval [— y / 2E/m, y/2E /m] with the restric- 
tion on the sum of the squares given by Eq. (|4j . In other 
words, the rescaled positions q with qi a = r.i a /L a are 
confined to the unit hypercube in dN dimensions, while 
the rescaled velocity components u = ym/ (2E) v are 
confined to the surface of the unit hypersphere in dN 
dimensions defined by the constraint u • u = 1. 

The state of the system is specified by the phase space 
vector of all velocities and positions T = (v,r), i.e., by 
2dN variables: the velocity components Vi a and the po- 
sition components ri a . However, these variables are not 
independent because of constraints. For spheres with 
random velocities and positions confined in a container 
with hard reflecting walls, the total energy E is conserved 
(microcanonical ensemble, constant NVE) and thus the 
number of independent variables is g — 2dN — 1. With 
periodic boundary conditions also the total linear mo- 
mentum 

N N 

P = ^2 m i v i ^m^Vi (5) 

i=l i=l 

and center of mass 

Ei=i m ^ l=1 

are conserved (molecular dynamics ensemble, constant 
NVEPG), and thus the number of independent variables 
drops to g = 2d{N - 1) - 1 = 2dN - 2d - 1. Symmetries 
in the positions and velocities may reduce g too; e.g. if all 
components i of T are pairwise symmetric with rispect 
to the origin, with both kinds of boundary conditions 
this point symmetry will stay on forever and g = dN — 1 
or g = 2d(N/2 - 1) - 1 = dN - 2d - 1 respectively. 
For the sake of simplicity, in presenting the theory we 
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will treat explicitly only the microcanonical case without 
symmetries. 

Following Khinchin, one can assume the uniform dis- 
tribution in the accessible portion of phase space as the 
starting point of statistical mechanics [24|. In our case, 
the measure of the accessible region of phase space is the 
product of the volume of the d./V-dimensional hypercube 
times the surface of the dTV-dimensional hypersphere, 



n = v 



2n dN/2 

T{dN/2) 



2E\ 
m J 



(dN-l)/2 



(7) 



This applies to a phase space whose coordinates are ve- 
locities and positions; of course the expression will be 
slightly different using momenta rather than velocities, 
or if the so-called density of states with respect to the 
energy, il' = dfl/dE, is used instead [24j . 

Khinchin's Ansatz is that the probability density func- 
tion (PDF) for points (v, r) in the permitted region of 
phase space is uniform. However, so far this has not 
been rigorously proved in general. In our case this Ansatz 
leads to the joint PDF for velocities and positions 

/v,r(x,y) = ^l{ x . x= 2£}(x)l n ^j_i^ ^j«(y), (8) 
where 1a(x) is the indicator function of the set A, 



1a(x) 



del' 



' \ if x e A 
if x i A 



(9) 



As the energy does not depend on positions, one can inte- 
grate over the latter, yielding a uniform PDF for particle 
velocities on the surface of a hypersphere, 



„ , . T(dN/2) ( m\(dN-i)/2 



(10) 



The marginalization of this joint PDF leads to the distri- 
butions of individual particle energies as well as of veloc- 
ity moduli and velocity components. To this purpose, it 
is convenient to study the relationship between Eq. (fTU)) 
and the symmetric Dirichlet distribution with parameter 
a. 

The PDF of the n-dimensional Dirichlet distribution 
with parameter vector a is 

/£(x;a) d = f 

^ n n 

n *r 1 n . 4=1 }(*). in) 

V ' i=l i=l 

Its value is zero outside the unit simplex 

S= jxe K" : Vi x l > OAj^x, = l| , (12) 
and thus Eq. pT|) can also be written 

n 

^"^Bssn^w- as) 

' ' i=l 



The normalization factor is given by the multinomial beta 
function, which can be defined through the gamma func- 
tion, 



B(a) = 



r(E?=i«0" 



(14) 



In the symmetric Dirichlet distribution all elements of 
the parameter vector a have the same value aj = a, 



F(a)] 



n 



l ls(x). 



(15) 



Notice that a = 1 gives the uniform distribution on S. 

It is convenient to work with adimensional variables. 
With the rescaling u ia = y / m/(2E)vi a introduced 
above, one gets the PDF 



, . , T(dN/2) , . 

/"W = 2lI dN/2 l{x-x=l}(x). 



A second transformation 



(16) 



(17) 



leads to a set of dN random variables each one with sup- 
port in [0,1] and such that 



N d 
i—l a—1 

The Jacobian for this transformation is 

JV d 



du 



9w 



odN n n 



-1/2 



(18) 



(19) 



i—l a—1 

Multiplying it by a factor 2 dN because each ±Ui a results 
in the same Wi a and by another factor 2 because of the 
constraint given by Eq. (|T5)) (for details see Song and 
Gupta Hi|), and replacing ^/n = T(l/2), the joint PDF 
of the variables Wi a can be expressed through the sym- 
metric Dirichlet PDF with parameter a = 1/2, 



/ w (x) = /°(x;l/2) 
T(dN/2) 



nn 



-1/2 



[r(i/2)]^iiii 

Now the normalized energy per particle, 



l s (x). (20) 



St = 



Ei 
E 



2E 



a=l 



(21) 



is the sum of d variables following the distribution given 
by Eq. ([20|) . As a consequence of the aggregation law for 
Dirichlet distributions, one finds that the joint PDF of 
all is 



/,(x) = / £ D (x;d/2) 
T{dN/2) 



[T(d/2)}" 



N 

n 



d/2-l 

V 1 



ls(x). 



(22) 
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It is interesting to notice that this is a uniform distribu- 
tion for d = 2; because of this, Boltzmann's 1868 method 
works in d = 2 dimensions, but fails in d — 3 dimensions 

The PDF of the normalized energies of single parti- 
cles can be obtained by a further marginalization of the 
symmetric Dirichlet distribution given by Eq. (|22p , using 
again the aggregation law. The result is a beta distribu- 
tion, whose PDF is 



f x (x;a,b) = 



1 



B(a,b) 



x^il-xf-H^ix); (23) 



recalling Eq. (HU, B(a,b) = T (a)T (b) /T (a + b) . Our case 
has the exponents a = d/2 and b = d(N — l)/2, 



d d(N — 1) 



x: —, 

' 2' 



T(dN/2) 



d/2-1 



T(d/2)T(d(N -l)/2)' 



(24) 



The transformation of variables Ei = Eei immediately 
leads to the PDF of particle energies, that follow a beta- 
Stacy distribution, 



x _ d d(N 
E' 2' 2 
r(dJV/2) 



1) 



d x 
6x E 

d/2 



T(d/2)T(d(N- l)/2 

jj; V d(JV-l)/2-l i 



(!) 



l[o, B ]W (25) 



for N > 1, and = 5(a: - £) for iV = 1. This result 

was obtained with a different method, without invoking 
the Dirichlet and beta distributions, by Shirts et al. [26l . 
Eq. (9)]. 

In the thermodynamic limit (N, V, E — > oo with 
N/V = p = constant and E/N = E = constant), 
Eq. (j25j) converges to a gam ma distribution, as discussed 
by Garibaldi and Scalas [23j, pages 121-122]. The gamma 
PDF is 



fx(x;a,b) 



del X" 



exp ( -- j l[ 0)OO )(ar). 



(26) 



b a T(a) 

A scale parameter is usually included in the definition of 
the gamma distribution, but it could always be set to 1 
absorbing it into the argument, 

f x ix;a,b) = fx (f ;M) \ = fx (f ;o) J- (27) 

Coming back to the thermodynamic limit of Eq. (1251) an- 
ticipated above, this is a gamma distribution with shape 
parameter a = d/2 and scale parameter b = 2E/(dm), 



d 2E\ 

dm\ d/2 x d / 2 - 1 

W 6XP 



/ dm V 
\2e) 



dmx\ 
^2EJ 



l[ 0jO o)(a!), (28) 



which is the familiar Boltzmann or Boltzmann-Gibbs dis- 
tribution for d = 2. 

The PDF of the velocity moduli, or speeds, of in- 
dividual particles can be obtained from fs 4 (x) replac- 
ing Vi = y/2Ei/m. The result is a transformed beta- 
Stacy distribution with the same exponents a = d/2 
and b — d(N — l)/2 as for the energies, but argument 
mx 2 /(2E), 



fvM = ft, 



p fmx 2 d d(N - 1) 



2E ' 2' 
T(dN/2) 



mx 



T(d/2)T(d(N- l)/2) 

d(N-l) . 



dx 2E ' 

d i 



2# 



x 1- 



2£ 



mi 



I (s) (29) 



for iV > 1, and f n (x) = 5{x - y / 2E/m) for N = 1. 
Also this result was obtained with a different method by 
Shirts et al. Hi. 



In the thermodynamic limit, Eq. (|29l) converges to the 
transformed gamma distribution with argument x 2 /2, 
shape parameter a = d/2 and scale parameter b = 
2E/(dm), 



, , . , , ^ d 2E\ d_x^ 
JvtW-Jv, 1 2 dm ) dx 2 



dm 
IF 



(x/2) d ~ 1 
W 



exp 



dmx' 
4E 



1 [o,oo)(a;) (30) 



which is the familiar Maxwell distribution for d = 3. 

The transformation from hyperspherical coordi- 
nates to cartesian coordinates v 2 — Yla=i v ia an< ^ 
(27r d / 2 /r(d/2))w^ 1 dw i = U d a =i dv ia leads from Eq. (J29]) 
to the PDF f Vi (x) of the single-particle velocity vectors 



/v,(x) 



T(dN/2) 



m \ 

T(d(N- l)/2) V2^bJ 

TO vd(iV-l)/2-" 
1 X • X 

2£ 



d/2 



L {x-x=2f} 



(x), (31) 



an equation obtained before too [26J, |2 

The direct marginalization [25J of the joint PDF of all 
velocities, Eq. (fT0|) . leads to the PDF f Via (x) of veloc- 
ity components, a result obtained integrating over all i 
except one and over all a except one. This is the quan- 
tity discussed by Maxwell [l|, and its derivation for any 
N is one of the main results in this paper. It turns out 
that the PDF of the velocity components is a transformed 
beta distribution with argument (l + ym/(2~B) x) /2 and 
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equal exponents a = b = (dN — l)/2, 

1 r PTE . /2El (^) 



1 / to 



2 \ 2E B((dN - l)/2, (dN - l)/2) 

1 I m x\ fl I m x 

2 + \l2E2J \2 ~ \I2E2 

i / r(div - i) 



(dN-3)/2 



2 dN-2 y 2ET 2 ((dN- l)/2) 



2E /2E 



2£ 



. ; ,1 /TrTx diV - 1 dN - 1 
/««„ I 2 + V2£2' 2 ' 2 



1 

dx V 2 



TO CC 

2~E 2 



(32) 



In the thermodynamic limit Eq. (|32[) converges to a nor- 
mal law with average \i = and variance c 2 = dE/(2m), 
i.e. the familiar Maxwell-Boltzmann distribution 



dnE 



exp 



(33) 



This is again related to a gamma distribution, since the 
positive half of the normal distribution can be expressed 



2ir<r 



exp 



— )l [0tOo) (x) = fZ[—;-,a )——. 



2 ' 2' 



da; 2 



In summary, all the known results for the relevant dis- 
tributions of the NVE ensemble can be obtained observ- 
ing that the normalized individual particle energies Si — 
Ei / E follow a symmetric multivariate Dirichlet distribu- 
tion with parameter a = 1/2 given by Eq. (|2T)]) . This is 
a direct consequence of the uniform-distribution assump- 
tion in Eq. ([8]) via a simple change of variables. Only for 
the velocity components it is necessary to marginalize the 
uniform distribution directly on the surface of the hyper- 
sphere and not on the simplex. Maxwell's Ansatz is vin- 
dicated by the fact that, in the thermodynamic limit, a 
normal distribution for velocity components is recovered, 
as well as their independence. Finally, for the NVEPG 
ensemble, the constraint given in Eq. ([5]) leads to differ- 
ent distributions for the relevant quantities introduced 
above. This will become clearer in the following. 



III. MOLECULAR DYNAMICS SIMULATIONS 

In molecular dynamics (MD) with continuous poten- 
tials, the equations of motion are integrated numerically 
using a constant time step; this approach is called time- 
driven. The larger the forces, the smaller the time step 
necessary to ensure energy conservation. With step po- 
tentials there are no forces acting on a distance, only 
impulsive ones at the exact time of impact. Therefore an 



event-driven approach is more appropriate: rather than 
until a fixed time step, the system is propagated until 
either the next collision or the next boundary crossing 

[IlH. 

The collision time tij between two particles i, j can be 
calculated from the mutual distance ry 



relative velocity Vj 



If h 



= r, 
vn ■ r 



Yj and the 
> the 



particles are moving away from each other and will not 
collide. Otherwise impact may happen at time tij when 
their distance becomes equal to the sum of their radii, 
i.e., ||ry + tyVyH = a. This is a second order problem 
with solutions 



(35) 



If the solutions are complex, no collision occurs. If the 
solutions are real, the smaller one, t~, corresponds to 
when the particles first meet, while the larger one, ij , to 
when they leave each other assuming they are allowed to 
interpenetrate. A negative collision time means that the 
event took place in the past. Because of the condition 
bij < 0, at least tfj > 0. If t~ < the particles overlap, 
which indicates an error. So the collision time is given 
by t~j 7 provided it is a positive real number. 

For a system of N hard spheres, at impact, assuming 
an elastic collision, the total kinetic energy E and the 
total linear momentum P are conserved (usually one sets 
P = at the beginning of the simulation by subtracting 
Y^Li Vj from each Vj). Assuming smooth surfaces, the 
impulse acts along the line of centers of the collision part- 
ners i and j given by fy; with equal masses, Vj changes 
to Vj + Avj and Vj changes to Vj — Av^ with 



bijTij 



1 y ) 1 y 



(36) 



where and Vy are evaluated at the instant of collision, 
and thus \\rij\\ = o. 

When a particle reaches a side of the unit box, peri- 
odic boundary conditions may require to "rebox" it by 
reintroducing it on the other side, while hard reflecting 
walls require to invert the velocity component perpen- 
dicular to the wall. After an event, be it a collision with 
another particle, a boundary crossing or a reflection at 
a boundary, the event calendar must be re-evaluatcd for 
pairs involving one of the event participants or a particle 
scheduled to collide with one of the event participants. 
All other particles are not influenced. Thus not every 
scheduled event actually takes place, because it can be 
invalidated by another earlier event, in which case it is 
erased from the priority queue. The latter is most com- 
monly handled by means of a binary tree (33] , which we 
realized with a multimap of the CH — V Standard Template 
Library [34j . The efficiency of this and alternative data 
structures for event scheduling has been analyzed exten- 
sively m, Hsj. 

The computational effort to search for min^j tij grows 
as the square of the number of particles. For large sys- 
tems it is advisable to divide the simulation box into 
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cells [37J. Provided cell boundary crossings are consid- 
ered too in the event list, two particles can collide only if 
they are located in the same cell or in adjacent cells. We 
chose cells with a side larger than a particle diameter. 
For more details on this and other algorithmic aspects in 
event-driven MD see Refs. [38l - l40l|. F or a parallel imple- 
mentation see Miller and Luding |4l| . 



V. NUMERICAL RESULTS 



IV. MONTE CARLO SIMULATIONS 



Except for especially ordered initial conditions, in- 
terparticle collisions computed by MD as explained in 
Sec. IIIII have mutual distance versors at collision fy uni- 
formly distributed on a unit half sphere in d dimensions 
such that, given relative velocities vy, the scalar prod- 
uct Vy ■ Tij is negative. Therefore the same distributions 
of velocities, and thus of derived quantities like energies, 
as in MD with periodic boundaries can be obtained by 
Monte Carlo (MC): after initializing the velocities of all 
hard spheres, the MC cycles consist in selecting a pair 
ij and a random versor fy such that Vy ■ fy < 0, and 



then in updating the velocities according to Eq. 
Hard reflecting walls can be included in the MC scheme 
by selecting with a certain frequency a sphere i and in- 
verting one of its velocity components Vi a . Altogether 
this is very much easier to code and faster to run, espe- 
cially for large numbers of particles N, than with MD, 
because no event list management is necessary. Moreover 
this scheme gives a useful insight into the mechanism of 
energy and momentum transfer. 

For a given initial state (a set of particle velocities), 
the MC dynamics denned above provides the realiza- 
tion of a Markov chain with symmetric transition ker- 
nel, meaning that P(v'|v) = P(v|v'), where v is the 
old velocity vector before the transition and v' is the 
new velocity vector after the transition. This Markov 
chain is homogeneous, as the transition probability does 
not depend on the time step. Invoking detailed balance, 
P(v'|v)P(v) = P(v|v')P(v'), the symmetry of the tran- 
sition kernel implies that the stationary distribution of 
this chain is uniform over the set of accessible states. If 
this set coincides with the surface of the velocity hyper- 
sphere, then the Markov chain is ergodic and one can 
hope to prove that the uniform distribution over the 
hypershere is also the equilibrium distribution for the 
Markov chain; see Sigurgeirsson [4^, Chapter 5] for the 
discussion of a related problem, and Meyn and Tweedie 
(43| for general methods. The results of MC simulations 
described below corroborate this conjecture and the al- 
gorithm outlined above is indeed an effective way of sam- 
pling the uniform distribution on the surface of a hyper- 
sphere. 



In Sec. HV} we have presented a simple stochastic model 
able to reproduce the same empirical equilibrium distri- 
bution for the random variables Vi a as obtained by MD 
simulations. In this sense, in principle, the above stochas- 
tic model and the related MC simulations can be used to 
estimate f Via {x) = lim^oo f Via (x, t). It turns out that 
the PDF f Vfa (x) for d = 2 and N = 2 is fitted by the 
arcsine law, 



1 



(37) 



ny/2E/m-x 2 
The name is due to its cumulative distribution function, 



x + 



1 



(38) 



If a Kolmogorov-Smirnov (KS) goodness-of-fit test (441 - 
I4H ] is performed comparing Eq. with the empirical 
cumulative distribution function of MC velocities for N — 
2, see Tab.Hl one gets that the null hypothesis of arcsine- 
distributed data cannot be rejected at the 5% significance 
level. For our sample size the value of the KS statistic 
is 3.2 x 1CT 4 with a p value of 0.99, the critical value 
being 9.6 x 10~ 4 . This test has not been repeated for 
MD velocities because the two empirical PDFs coincide 
as shown in Fig. [5] 



N 



Kolmogorov-Smirnov 



Pks 



2 

3 

10 
100 
1000 
10 000 



4.1 x 10" 

1.0 x 10" 
1.3 x 10" 

7.1 x 10" 

6.2 x 10" 
1.1 x 10" 



0.99 
0.24 
0.07 
0.69 
0.84 
0.15 



TABLE I: Comparison between the empirical probability den- 
sity functions of the velocity components from MC with pe- 
riodic boundary conditions and d = 2 by means of the KS 
test. The sample size is 2 x 10 6 for every value of N. At the 
5% significance level the critical value is 9.6 x 10 -4 . The null 
hypothesis of equally distributed data can never be rejected. 

For d — 2 and N — 3, f Via (x) is described by the 
semicircle law, 



/««, 0*0 



m 
2^E 



AE 



(39) 



Its cumulative distribution function again contains an 
arcsine, 



mx / m 

F v . (x) = =■* — =■ — x z 

1Q V ; AttE V AE 



1 / Pm \ 1 

(?o) 

If a KS goodness-of-fit test is performed comparing 
Eq. (|40p with the empirical cumulative distribution func- 
tion of MC velocities for N = 3, one gets that the null 
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FIG. 1: (Color online) Probability density functions of the velocity components (top), the velocity modulus (middle) and the 
energy (bottom) for d — 2 (left) and d = 3 (right) with E — 1, hard reflecting walls and zero total momentum. Lines: theory; 
empty symbols: MD; full symbols: MC. Delta functions are made visible by a vertical line for the theory and by rescaling down 
to 1 the data point that would otherwise be out of scale. 



hypothesis of semicircle-distributed data cannot be re- 
jected at the 5% significance level. For our sample size 
the value of the KS statistic is 6.0 x 10~ 4 with a p value 
of 0.46, the critical value being 9.6 x 10 -4 . This test 
has not been repeated for MD velocities because the two 
empirical PDFs coincide as shown in Fig. [2l 

For d = 3 and N = 2, f Via (x) is a uniform distribution 

on[-AV2]- 

The arcsine distribution (d = 2, N — 2) is given by 
Eq. (|32p with a = 1/2; the semicircle distribution (d = 



2, N = 3) is given by a — 3/2; the case with d = 2, N — 
4 is given by a = 5/2; the uniform distribution (d — 

3, N = 2) is given by a = 1; the case with d = 3, N = 3 
is given by a = 5/2; and so on. 

The empirical density f Via (x) is well approximated by 
a normal law already for N = 1000 hard disks, as shown 
in Tab. [TT1 where the results of two non-parametric tests 
for normality, Lilliefors 47J and Jarque and Bera [48|, |49| , 
are presented when N = 10, 100, 1000, 10 000; again these 
tests are done only for MC velocities. 
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FIG. 2: (Color online) Probability density functions of the velocity components (top), the velocity modulus (middle) and the 
energy (bottom) for d — 2 (left) and d = 3 (right) with E = 1, periodic boundary conditions and zero total momentum. Lines: 
theory; empty symbols: MD; full symbols: MC. Notice that the velocity components and energy MD data for N = 3 deviate 
slightly from the theory and the MC data. Delta functions are made visible by a vertical line for the theory and by rescaling 
down to 1 the data point that would otherwise be out of scale. 



VI. DISCUSSION AND CONCLUSIONS 

To summarize what we have done, in a system of N 
hard balls in a d-dimensional volume V the velocity com- 
ponents, the velocity modulus and the energies of the 
spheres or disks are well reproduced by transformed beta 
distributions with different arguments and shape param- 
eters depending on N, d, the total energy E, and the 
boundary conditions; in the thermodynamic limit these 



distributions converge to transformed gamma distribu- 
tions with different arguments and shape or scale pa- 
rameters, corresponding respectively to the Gaussian, i.e. 
Maxwell-Boltzmann, the Maxwell, and the Boltzmann 
or Boltzmann-Gibbs distribution. We showed this the- 
oretically using Khinchin's Ansatz, and performed sta- 
tistical goodness-of-fit tests on systematic MD and MC 
computer simulations of an increasing number N of hard 
disks or spheres starting from 2 in the microcanonical 
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1 V 


1 ill f ^ f /" ^ T"0 


PL 


J ell LJUC- JDCI a, 


PJB 


10 


0.008* 


< io- 3 


7.4 x 10 J * 


< 10 _J 


100 


7.36 x 10" 4 * 


0.01 


29.6* 


< 10~ 3 


1000 


4.79 x 10" 4 


0.35 


5.70 


0.06 


10 000 


4.37 x 10 -4 


0.50 


2.36 


0.31 



TABLE II: Results of two non-parametric normality tests for 
the empirical probability density function of the velocity com- 
ponents from MC with periodic boundary conditions when 
d — 2: Lilliefors (L) and Jarque and Bera (JB). The sample 
size is n = 2 x 10 6 . At the 5% significance level the criti- 
cal value is 6.43 x 10~ 4 for the L test and 5.99 for the JB 
test. The star indicates that the null hypothesis of normally 
distributed data can be rejected. 

ensemble (constant NVE, hard reflecting walls) and in 
the molecular dynamics ensemble (constant NVEPG, 
periodic boundary conditions). Even if these results are 
not entirely new, we presented comprehensively both the 
analytical derivations and the numerical checks, the lat- 
ter both by MD and MC; moreover, we realized that all 
these distributions are actually variants of the beta or the 
gamma distribution, and can be derived from the Dirich- 
let distribution; and we pointed out that for values of N 
as low as 2 or 3 the shapes of these distributions can be 
quite different from those in the thermodynamic limit; in 
particular, they can become uniform or even bimodal. 

The MD simulations presented above corroborate 
Boltzmann's ergodic hypothesis [50| both for the NVE 
and the NVEPG ensembles. The proofs of ergodicity for 
similar systems used the so-called Chernov-Sinai Ansatz, 
namely the almost sure hyperbolicity of singular orbits 
[5l| ; hyperbolicity means a non-zero Lyapunov exponent 
almost everywhere with respect to the Liouville measure. 
It was Sinai who, earlier [52J, had updated Boltzmann's 
ergodic hypothesis. One should prove that every hard- 
ball system on a flat torus is fully hyperbolic and er- 
godic, after fixing its total energy, momentum, and cen- 
ter of mass. This rephrasing of Boltzmann's hypothe- 
sis is known as the Boltzmann-Sinai ergodic hypothesis. 



More recently, Simanyi proved the Boltzmann-Sinai er- 
godic hypothesis in full generality for hard ball systems 

MD simulations show that Khinchin's Ansatz is justi- 
fied for systems of hard balls. We would like to stress 
that this is a consequence of the microscopic dynamics 
and not of any a priori maximum-entropy principle. The 
uniform distribution on the accessible phase-space region 
is indeed the maximum-entropy distribution. Therefore, 
maximum-entropy methods do work well and all the dis- 
tributions in Sec. [TT] could be obtained by maximum- 
entropy methods: the beta and gamma distributions are 
actually the maximum-entropy distributions with given 
first moment, and possibly some other constraint, on a 
finite and a semi-infinite interval respectively. However, 
this is so only because the dynamics uniformly samples 
the accessible phase-space region and not the other way 
round. In different frameworks, e.g. in biology or eco- 
nomics, maximum-entropy assumptions might lead to 
wrong results for the equilibrium distribution of a sys- 
tem, if its dynamics is not specified or carefully studied. 

The distributions derived in Sec. [TT] are a benchmark 
for random partition models popular in econophysics. 
Pure exchange models often lead to the same distribu- 
tions mum. 
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